tau=.1;    % time step

ggap=70;

T=100;
time=tau:tau:T;

type=[];
type(end+1)=-1;      % soma
type(end+(1:2))=0;   % trunk of apical dendrite 
type(end+(1:2))=1;   % C_1
type(end+(1:6))=2;   % C_2
type(end+(1:12))=3;  % C_3
type(end+(1:24))=4;  % C_4
n=length(type);      % the number of compartments

m=max(type)-min(type)+1;  % the number of compartment types

mvr=-60+zeros(1,m);
mvt=[-50,-50+zeros(1,m-1)];
mp=[3,3+zeros(1,m-1)];
mC=[100,100+zeros(1,m-1)];
ma=[0.01,0.01+zeros(1,m-1)];
mb=[5,5+zeros(1,m-1)];
mc=[-55,-35+zeros(1,m-1)];
md=[500,1000+zeros(1,m-1)];
mspike_peak=[50,10+zeros(1,m-1)];

vr=mvr(2+type);
vt=mvt(2+type);
p=mp(2+type);
C=mC(2+type);
a=ma(2+type);
b=mb(2+type);
c=mc(2+type);
d=md(2+type);

spike_peak=mspike_peak(2+type);

mgapm=ggap*[1 1 1 1 1 1];      % conductance (sensitivity) to j-1 (mother) compartment
mgapd=ggap*[1 1 1 1 1 1];     % conductance (sensitivity) to j+1 (daughter) compartment
gapm=mgapm(2+type);
gapd=mgapd(2+type);


nbrs=zeros(n,3);
nbrs(1,:)=[1 1 2];     % self-gap junctions (no effect), C_0
nbrs(2,:)=[1 2 3];     % C_s, self (no effect), C_0
nbrs(3,:)=[2 4 5];     % C_0, C_1, C_1
nbrs(4,:)=[3 6 7];     % C_0, C_2, C_2
nbrs(5,:)=[3 8 9];     % C_0, C_2, C_2
nbrs(6,:)=[4 12 13];   % C_1, C_3, C_3
nbrs(7,:)=[4 14 15];   % C_1, C_3, C_3
nbrs(8,:)=[5 16 17];   % C_1, C_3, C_3
nbrs(9,:)=[5 18 19];   % C_1, C_3, C_3
nbrs(10,:)=[1 20 21];   % C_s, C_3, C_3
nbrs(11,:)=[1 22 23];   % C_s, C_3, C_3
nbrs(12,:)=[6 24 25];   % C_2, C_4, C_4
nbrs(13,:)=[6 26 27];   % C_2, C_4, C_4
nbrs(14,:)=[7 28 29];   % C_2, C_4, C_4
nbrs(15,:)=[7 30 31];   % C_2, C_4, C_4
nbrs(16,:)=[8 32 33];   % C_2, C_4, C_4
nbrs(17,:)=[8 34 35];   % C_2, C_4, C_4
nbrs(18,:)=[9 36 37];   % C_2, C_4, C_4
nbrs(19,:)=[9 38 39];   % C_2, C_4, C_4
nbrs(20,:)=[10 40 41];   % C_2, C_4, C_4
nbrs(21,:)=[10 42 43];   % C_2, C_4, C_4
nbrs(22,:)=[11 44 45];   % C_2, C_4, C_4
nbrs(23,:)=[11 46 47];   % C_2, C_4, C_4
nbrs(24,:)=[12 24 24];   % C_3, self, self
nbrs(25,:)=[12 25 25];   % C_3, self, self
nbrs(26,:)=[13 26 26];   % C_3, self, self
nbrs(27,:)=[13 27 27];   % C_3, self, self
nbrs(28,:)=[14 28 28];   % C_3, self, self
nbrs(29,:)=[14 29 29];   % C_3, self, self
nbrs(30,:)=[15 30 30];   % C_3, self, self
nbrs(31,:)=[15 31 31];   % C_3, self, self
nbrs(32,:)=[16 32 32];   % C_3, self, self
nbrs(33,:)=[16 33 33];   % C_3, self, self
nbrs(34,:)=[17 34 34];   % C_3, self, self
nbrs(35,:)=[17 35 35];   % C_3, self, self
nbrs(36,:)=[18 36 36];   % C_3, self, self
nbrs(37,:)=[18 37 37];   % C_3, self, self
nbrs(38,:)=[19 38 38];   % C_3, self, self
nbrs(39,:)=[19 39 39];   % C_3, self, self
nbrs(40,:)=[20 40 40];   % C_3, self, self
nbrs(41,:)=[20 41 41];   % C_3, self, self
nbrs(42,:)=[21 42 42];   % C_3, self, self
nbrs(43,:)=[21 43 43];   % C_3, self, self
nbrs(44,:)=[22 44 44];   % C_3, self, self
nbrs(45,:)=[22 45 45];   % C_3, self, self
nbrs(46,:)=[23 46 46];   % C_3, self, self
nbrs(47,:)=[23 47 47];   % C_3, self, self


v=zeros(length(time),n);
v(1,:)=vr;
u=0*v;
g_ampa=0+0*v;
I=0+0*v;
st=65/tau;

decay = exp(-tau/5);

compartments = [2 3 5 9];
   g_ampa(floor(1/tau), 2 )=6.0;
   g_ampa(floor(30/tau), 3 )=6.0;
   g_ampa(floor(60/tau), 5 )=6.0;
   g_ampa(floor(90/tau), 9 )=6.0;
   
%   C(1:17) = [100 85 70 55  55 40 40 40 40 25 25 25 25 25 25 25 25 ] ;
%   p(1:17) = [3   2.5 2 1.5 1.5 1 1  1  1  1   1  1  1  1  1  1  1 ] ;
   
   
   for i=2:length(time)
       g_ampa(i,compartments)=g_ampa(i-1,compartments)*decay+g_ampa(i,compartments);
%       g_ampa(i,25)=300/2*t*exp(-2*t);
   end;

for i=1:length(time)-1    
 dv=p.*(v(i,:)-vr).*(v(i,:)-vt) -u(i,:) + I(i,:) + g_ampa(i,:).*(0-v(i,:)) + ...
     gapm.*(v(i,nbrs(:,1))-v(i,:)) + ...
     gapd.*(v(i,nbrs(:,2))+v(i,nbrs(:,3))-2*v(i,:));
 v(i+1,:)=v(i,:)+tau*dv./C;
 u(i+1,:)=u(i,:)+tau*(a.*(b.*(v(i,:)-vr)-u(i,:)));
 fired =find(v(i+1,:)>spike_peak);
 v(i,fired)=spike_peak(fired);
 v(i+1,fired)=c(fired);
 u(i+1,fired)=u(i+1,fired)+d(fired);
end;

plot(v(:,compartments),'.')

legend(num2str(compartments(1)),num2str(compartments(2)),num2str(compartments(3)),num2str(compartments(4)))
